# example 2: compare mean and median of normal distribution

n <- 100
nsim <- 10000

xmean=NULL
xmedian=NULL

for(i in 1:nsim){
x=rnorm(n,1,1)
xmean[i]=mean(x)
xmedian[i]=median(x)
}

mse_mean=mean((xmean-1)**2)
mse_median=mean((xmedian-1)**2)
cat('mse of mean=', mse_mean, sep="\n", 'mse of median=', mse_median,fill=T)

cat('mean of mean=',mean(xmean), sep="\n", 'mean of median=', mean(xmedian),fill=T)


hist(xmean)
hist(xmedian)

# lets include a few outliers
# after x=rnorm(n,1,1) command include the following (that is the laziest 
# way of creating outliers, for more professional work think about the 
# composition method we covered)

x[1:5]=10

# Example 4:

n=100
nsim=1000

mu=seq(-1,3,0.05)
mu

for(j in 1:length(mu)){

total=0

for(i in 1:nsim){
x=rnorm(n,mu[j],1)
xmean=mean(x)
if(0.8<xmean&xmean<1.2) total=total+1
}
power_test=1-total/nsim
write.table(cbind(mu[j],power_test),"sim.ex.4.csv",sep=",",row.names=F,col.names=F,append=T)
}

results=read.table("sim.ex.4.csv",sep=",")

head(results)
dim(results)

plot(results[,1],results[,2],xlab="mu",ylab="Prob of rejecting Ho",type="n")
lines(results[,1],results[,2])

# estimated alpha:
results[results[,1]==1,]